#include <iostream>
#include <vector>
#include <cmath>
#include <string>
#include <fstream>
#include <sstream> // 用于动态生成文件名
#include "CubicBSplines.h"
#define PI 3.14159265358979323846

// 定义函数类型
typedef std::vector<double> (*ParametricFunction)(double);

// 定义原始参数方程 f1
std::vector<double> f1(double t) {
    std::vector<double> result(2);
    result[0] = sqrt(3.0) * sin(2 * PI * t);
    result[1] = 2.0 * (sqrt(3.0) * cos(2 * PI * t) + sqrt(abs(sqrt(3.0) * sin(2 * PI * t)))) / 3.0;
    return result;
}

// 定义 f2
std::vector<double> f2(double t) {
    std::vector<double> result(2);
    result[0] = sin(t) + t * cos(t);
    result[1] = cos(t) - t * sin(t);
    return result;
}

// 定义 f3 (三维曲线)
std::vector<double> f3(double t) {
    std::vector<double> result(3);
    result[0] = sin(cos(t)) * cos(sin(t));
    result[1] = sin(cos(t)) * sin(sin(t));
    result[2] = cos(cos(t));
    return result;
}
std::pair<double, double> generateBoundaryConditions(const std::vector<double>& knots, const std::vector<double>& values) {
    size_t n = knots.size();
    if (n < 3) {
        throw std::invalid_argument("Not enough points for higher-order differences. At least 3 points are required.");
    }

    // 计算 m_1 的高阶差分
    double h1 = knots[1] - knots[0];
    double h2 = knots[2] - knots[1];
    double d1 = (values[1] - values[0]) / h1;   // 一阶差分
    double d2 = (values[2] - values[1]) / h2;   // 一阶差分
    double m_1 = d1 + (d2 - d1) * h1 / (h1 + h2); // 插值二阶导数外推

    // 计算 m_n 的高阶差分
    double hn_1 = knots[n - 1] - knots[n - 2];
    double hn_2 = knots[n - 2] - knots[n - 3];
    double dn_1 = (values[n - 1] - values[n - 2]) / hn_1; // 一阶差分
    double dn_2 = (values[n - 2] - values[n - 3]) / hn_2; // 一阶差分
    double m_n = dn_1 + (dn_1 - dn_2) * hn_1 / (hn_1 + hn_2); // 插值二阶导数外推

    return {m_1, m_n};
}

// 处理三维拟合并绘制
void plot3DCurveFitting(const std::vector<double>& knots, double start, double end, ParametricFunction func,
                        const std::string& dataFilename, const std::string& scriptFilename, const std::string& outputFilename) {
    std::ofstream outFile(dataFilename);
    std::ofstream scriptFile(scriptFilename);

    std::vector<double> x_values, y_values, z_values;

    // 计算原始数据点
    for (double t : knots) {
        std::vector<double> point = func(t);
        x_values.push_back(point[0]);
        y_values.push_back(point[1]);
        z_values.push_back(point[2]);
    }
    std::pair<double, double> m_1_n = generateBoundaryConditions(knots, x_values);
    std::pair<double, double> m_1_n2 = generateBoundaryConditions(knots, y_values);
    std::pair<double, double> m_1_n3 = generateBoundaryConditions(knots, z_values);
    double m_1_x = m_1_n.first;
    double m_n_x = m_1_n.second;
    double m_1_y = m_1_n2.first;
    double m_n_y = m_1_n2.second;
    double m_1_z = m_1_n3.first;
    double m_n_z = m_1_n3.second;
    // 生成样条实例（假设已有 CubicBSpline_natural 类支持）
    CubicBSpline_natural spline_x(knots, x_values);
    CubicBSpline_natural spline_y(knots, y_values);
    CubicBSpline_natural spline_z(knots, z_values);
    CubicBSpline_complete spline_x2(knots, x_values,m_1_x,m_n_x);
    CubicBSpline_complete spline_y2(knots, y_values,m_1_y,m_n_y);
    CubicBSpline_complete spline_z2(knots, z_values,m_1_z,m_n_z);
    CubicBSpline_periodic spline_x3(knots, x_values);
    CubicBSpline_periodic spline_y3(knots, y_values);
    CubicBSpline_periodic spline_z3(knots, z_values);

    // 生成样条曲线数据
    for (int i = 0; i < 200; ++i) {
        double ti = start + (end - start) * i / 199.0;
        double xi = spline_x.B_value(ti);
        double yi = spline_y.B_value(ti);
        double zi = spline_z.B_value(ti);
        double xi2 = spline_x2.B_value(ti);
        double yi2 = spline_y2.B_value(ti);
        double zi2 = spline_z2.B_value(ti);
        double xi3 = spline_x3.B_value(ti);
        double yi3 = spline_y3.B_value(ti);
        double zi3 = spline_z3.B_value(ti);
        outFile << xi << " " << yi << " " << zi << " 2" << std::endl;
        outFile << xi2 << " " << yi2 << " " << zi2 << " 3" << std::endl;
        outFile << xi3 << " " << yi3 << " " << zi3 << " 4" << std::endl; // 第四列标识为样条曲线
    }

    // 绘制原始函数
    for (double t = start; t <= end; t += 0.01) { // 使用适当的步长
        std::vector<double> point = func(t);
        outFile << point[0] << " " << point[1] << " " << point[2] << " 1" << std::endl; // 原始函数数据
    }

    outFile.close();

    // 创建 gnuplot 脚本
    scriptFile << "set terminal pngcairo enhanced color size 1200,800\n";
    scriptFile << "set output '" << outputFilename << "'\n";
    scriptFile << "set title '3D Curve Fitting Using B-Splines'\n";
    scriptFile << "set xlabel 'x'\n";
    scriptFile << "set ylabel 'y'\n";
    scriptFile << "set zlabel 'z'\n";
    scriptFile << "set grid\n";
    scriptFile << "set key outside\n";
    scriptFile << "splot \\\n";
    scriptFile << "     '< grep \" 1$\" " << dataFilename << "' using 1:2:3 with lines title 'Original Function' lc rgb 'red' lw 2, \\\n";
    scriptFile << "     '< grep \" 2$\" " << dataFilename << "' using 1:2:3 with lines title 'Natural B-Spline' lc rgb 'blue' lw 2,\\\n";
    scriptFile << "     '< grep \" 3$\" " << dataFilename << "' using 1:2:3 with lines title 'Complete B-Spline' lc rgb 'green' lw 2,\\\n";
    scriptFile << "     '< grep \" 4$\" " << dataFilename << "' using 1:2:3 with lines title 'Periodic B-Spline' lc rgb 'orange' lw 2\n";
    scriptFile.close();

    // 执行 gnuplot 脚本
    int returnCode = system(("gnuplot " + scriptFilename).c_str());
    if (returnCode == 0) {
        std::cout << "3D Plot saved as '" << outputFilename << "'" << std::endl;
    } else {
        std::cerr << "Error: Failed to execute gnuplot script." << std::endl;
    }
}
// 生成均匀分布的节点
void generateKnots_uniform(std::vector<double>& knots, double start, double end, int n) {
    knots.resize(n);
    double step = (end - start) / (n - 1);
    for (int i = 0; i < n; i++) {
        knots[i] = start + i * step;
    }
}

// 计算两点间的欧几里得距离
double euclideanDistance(const std::vector<double>& p1, const std::vector<double>& p2) {
    double sum = 0.0;
    for (size_t i = 0; i < p1.size(); ++i) {
        sum += (p2[i] - p1[i]) * (p2[i] - p1[i]);
    }
    return std::sqrt(sum);
}

// 使用累积弦长（Cumulative Chordal Length）生成节点
void generateKnots_cumulative(std::vector<double>& knots, const std::vector<std::vector<double>>& points) {
    size_t n = points.size();
    knots.resize(n);

    // 累积弦长
    knots[0] = 0.0; // 第一个节点从 0 开始
    for (size_t i = 1; i < n; ++i) {
        knots[i] = knots[i - 1] + euclideanDistance(points[i - 1], points[i]);
    }

    // 将累积弦长归一化到 [0, 1] 范围
    double totalLength = knots.back();
    for (size_t i = 0; i < n; ++i) {
        knots[i] /= totalLength;
    }
}

// 绘制曲线拟合并存储到单独文件
void plotCurveFitting(const std::vector<double>& knots, double start, double end, ParametricFunction func,
                      const std::string& dataFilename, const std::string& scriptFilename, const std::string& outputFilename) {
    std::ofstream outFile(dataFilename);
    std::ofstream scriptFile(scriptFilename);

    std::vector<double> x_values, y_values;

    // 计算原始数据点
    for (double t : knots) {
        std::vector<double> point = func(t);
        x_values.push_back(point[0]);
        y_values.push_back(point[1]);
    }
std::pair<double, double> m_1_n = generateBoundaryConditions(knots, x_values);
    std::pair<double, double> m_1_n2 = generateBoundaryConditions(knots, y_values);
    double m_1_x = m_1_n.first;
    double m_n_x = m_1_n.second;
    double m_1_y = m_1_n2.first;
    double m_n_y = m_1_n2.second;
    // 生成样条实例（假设已有 CubicBSpline_natural 类支持）
    CubicBSpline_natural spline_x(knots, x_values);
    CubicBSpline_natural spline_y(knots, y_values);
    CubicBSpline_complete spline_x2(knots, x_values,m_1_x,m_n_x);
    CubicBSpline_complete spline_y2(knots, y_values,m_1_y,m_n_y);
    // 生成样条曲线数据
    for (int i = 0; i < 200; ++i) {
        double ti = start + (end - start) * i / 199.0;
        double xi = spline_x.B_value(ti);
        double yi = spline_y.B_value(ti);
        double xi2 = spline_x2.B_value(ti);
        double yi2 = spline_y2.B_value(ti);
        outFile << xi << " " << yi  << " 2" << std::endl;
        outFile << xi2 << " " << yi2  << " 3" << std::endl;
    }

    // 绘制原始函数
    for (double t = start; t <= end; t += 0.01) { // 使用适当的步长
        std::vector<double> point = func(t);
        outFile << point[0] << " " << point[1] << " " << point[2] << " 1" << std::endl; // 原始函数数据
    }

    outFile.close();

     // 创建 gnuplot 脚本
    scriptFile << "set terminal pngcairo enhanced color size 1200,800\n";
    scriptFile << "set output '" << outputFilename << "'\n";
    scriptFile << "set title '2D Curve Fitting Using B-Splines'\n";
    scriptFile << "set xlabel 'x'\n";
    scriptFile << "set ylabel 'y'\n";
    scriptFile << "set grid\n";
    scriptFile << "set key outside\n";
    scriptFile << "plot \\\n";
    scriptFile << "     '< grep \" 1$\" " << dataFilename << "' using 1:2 with lines title 'Original Function' lc rgb 'red' lw 2, \\\n";
    scriptFile << "     '< grep \" 2$\" " << dataFilename << "' using 1:2 with lines title 'Natural B-Spline' lc rgb 'blue' lw 2, \\\n";
    scriptFile << "     '< grep \" 3$\" " << dataFilename << "' using 1:2 with lines title 'Complete B-Spline' lc rgb 'green' lw 2,\\\n";
    scriptFile.close();

    // 执行 gnuplot 脚本
    int returnCode = system(("gnuplot " + scriptFilename).c_str());
    if (returnCode == 0) {
        std::cout << "Plot saved as '" << outputFilename << "'" << std::endl;
    } else {
        std::cerr << "Error: Failed to execute gnuplot script." << std::endl;
    }
}
int main() {
    std::vector<ParametricFunction> functions = {f1, f2, f3};
    std::vector<std::string> functionNames = {"f1", "f2", "f3"};
    std::vector<std::pair<double, double>> t_ranges = {{0, 1}, {0, 6 * PI}, {0, 2 * PI}};
    std::vector<int> N_values = {10, 40, 160}; // 节点数量

    for (size_t funcIdx = 0; funcIdx < functions.size(); ++funcIdx) {
        ParametricFunction func = functions[funcIdx];
        std::string funcName = functionNames[funcIdx];
        double t_start = t_ranges[funcIdx].first;
        double t_end = t_ranges[funcIdx].second;

        for (int N : N_values) {
            std::cout << "Processing " << funcName << " with N = " << N << "..." << std::endl;

            // 使用均匀分布生成节点
            std::vector<double> knots_uniform;
            generateKnots_uniform(knots_uniform, t_start, t_end, N);

            // 动态生成文件名（均匀分布节点）
            std::ostringstream uniform_data, uniform_script, uniform_output;
            uniform_data << funcName << "_uniform_data_N" << N << ".txt";
            uniform_script << funcName << "_uniform_plot_N" << N << ".gp";
            uniform_output << funcName << "_uniform_curve_N" << N << ".png";

            // 对 f3 使用 3D 绘图，否则绘制 2D
            if (funcName == "f3") {
                plot3DCurveFitting(knots_uniform, t_start, t_end, func, uniform_data.str(), uniform_script.str(), uniform_output.str());
            } else {
                plotCurveFitting(knots_uniform, t_start, t_end, func, uniform_data.str(), uniform_script.str(), uniform_output.str());
            }

            // 使用累积弦长生成节点
            std::vector<std::vector<double>> points(N);
            for (int i = 0; i < N; ++i) {
                double t = t_start + (t_end - t_start) * i / (N - 1);
                points[i] = func(t);
            }

            std::vector<double> knots_cumulative;
            generateKnots_cumulative(knots_cumulative, points);

            // 动态生成文件名（累积弦长节点）
            std::ostringstream cumulative_data, cumulative_script, cumulative_output;
            cumulative_data << funcName << "_cumulative_data_N" << N << ".txt";
            cumulative_script << funcName << "_cumulative_plot_N" << N << ".gp";
            cumulative_output << funcName << "_cumulative_curve_N" << N << ".png";

            // 对 f3 使用 3D 绘图，否则绘制 2D
            if (funcName == "f3") {
                plot3DCurveFitting(knots_cumulative, t_start, t_end, func, cumulative_data.str(), cumulative_script.str(), cumulative_output.str());
            } else {
                plotCurveFitting(knots_cumulative, t_start, t_end, func, cumulative_data.str(), cumulative_script.str(), cumulative_output.str());
            }
        }
    }

    return 0;
}